function [r] = speccalc2(dfile, v, v2, time1, time2)
% compute power spectrum of v data from t1 to t2
% from data first
% Based on Speccalc (1-2 Mar 1999. P. Manis)
% 5/3/99 P. Manis
% basic results of analysis are stored in the result structure
% as:
% result.filename = 'filename'
% result.records = 'frec - lrec'
% result.f0 = peak osc. freq from spectral analysis
% result.max = psd at f0
%
%
r = []; % result structure...
disp(sprintf('speccalc2 with: t1= %6.1f  t2=%6.1f', ...
   time1, time2))
t = make_time(dfile);
%[v, t] = datac('getv'); % get the data first
%dfile = datac('getdfile');
r.filename = dfile.filename;
r.records = sprintf('%d - %d', min(dfile.record), max(dfile.record));
rate = dfile.rate*dfile.nr_channel/1000;
fsamp = 1000/rate; % get sampling frequency
fco = 250;		% cutoff frequency in Hz
wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
if(wco < 1) % if wco is > 1 then this is not a filter!
   [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
   vsmo = DigitalFilt(b, a, v')'; % filter all the traces...
	vsmo2 = DigitalFilt(b, a, v2')';
else
   vsmo = v;
end
v = vsmo;
v2 = vsmo2;

t1=min(find(t>time1));
t2=min(find(t>time2));
nr = size(v);
nfft =2048;
h = findobj('Tag', 'SpecPlot'); % check for pre-existing window
if(isempty(h)) % if none, make one
   h = figure('Tag', 'SpecPlot', 'Name', ...
      sprintf('Spectral Analysis of Traces'), 'NumberTitle', 'off');
end
figure(h); % otherwise, select it
clf; % always clear the window...
%f = (1000/rate)*(0:(nfft/2))/nfft;
subplot(3,1,1);
title(sprintf('Spectral Analysis for %s Records %d-%d', dfile.filename, min(dfile.record), max(dfile.record)));
ylabel('mV');
xlabel('ms');
tm=make_time(dfile);
plot(tm, v);
Pv=zeros(nr(1), nfft);
ngrand = 0;
novlap = round(2/rate);
windlen = round(50/rate);

nharm = 24;
r.f60=60*[1:nharm]; % harmonics of 60 Hz.
for j = 1:nr(1) % compute power spectrum for each of the records individually
   %  v(j,:) = v(j,:) -mean(v(j,:));
   [Pv, f] = spect(v(j,t1:t2), nfft, novlap, windlen, 1000/rate, 'linear');
   r.f=f;
   [Pv2, f] = spect(v2(j,t1:t2), nfft, novlap, windlen, 1000/rate, 'linear');
   for i=1:length(r.f60)
      f0=r.f60(i);
      [q, fl1]=min(abs(f-f0)); % fl(1);
      r.a(i)=Pv(fl1); % these are the amplitudes at each of the harmonics
      r.b(i)=Pv2(fl1);
   end
   r.a0= Pv(1);
   r.b0= Pv2(1);
   r.a  = r.a / r.a0;
   r.b = r.b / r.b0;
   vsmo3(j,:) = vsmo2(j,:) - (mean(r.a(2))/mean(r.b(2)))*vsmo(j,:);
   subplot(3,1,2);
   loglog(f, sqrt(Pv(1:((nfft/2)+1))), '-k');
   hold on;
   loglog(r.f60, sqrt(r.a), 'xr');
   loglog(f, sqrt(Pv2(1:((nfft/2)+1))), '-b');
   axis([10 5000 0.1 3000]);
   r.a = r.a/fsamp;
   
   novlap = round(2/rate);
   windlen = round(50/rate);
   % now compute the cosine equivalent of the harmonic series... based on the average...
   % first we correct the amplitudes so that the correspond to a point frequency...
   ssin = 0*(t1:t2);
   for i=1:length(r.f60)
      f0=r.f60(i);
      ssin = ssin + r.a(i)*sin((2*pi*f0/1000)*t(t1:t2)); % get the angles at 0 frequency
   end
   X = fft(vsmo(j,t1:t2));
   XP = fft(ssin);
   ph = angle(X(1:length(X/2)))-angle(XP(1:length(XP/2)));
   Y=ifft(X);
   ssin=0*(t1:t2);
   for i = 1:length(r.f60)
      f0 = r.f60(i);
      [q, fl1] = min(abs(f-f0));
      r.phm(j,i) = mean(ph(fl1));
      ssin = ssin + r.a(i)*sin(((2*pi*f0/1000)*t(t1:t2))-r.phm(j,i)); % get mean phase in here too.
   end
   r.phm(1:max(length(r.phm), min(3, length(r.phm))))
   subplot(3,1,3)
 %  plot(t(t1:t2), ssin);
 plot(t(t1:t2), mean( vsmo2(:, (t1:t2))));
 newsig(j,:) = vsmo(j,t1:t2) - ssin;
   subplot(3,1,1);
    plot(t(t1:t2), vsmo3(j,(t1:t2)), '-r', t(t1:t2), vsmo2(j,(t1:t2)), '-k');
   pause
end
%subplot(3,1,1);
%plot(t(t1:t2), newsig);

return;

